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Abstract: We investigate the parametrization issue for discrete-time stable all- 
pass multivariable systems by means of a Schur algorithm involving a Nudelman 
interpolation condition. A recursive construction of balanced realizations is as- 
sociated with it, that possesses a very good numerical behavior. Several atlases 
of charts or families of local parametrizations are presented and for each atlas 
a chart selection strategy is proposed. The last one can be viewed as a nice 
mutual encoding property of lossless functions and turns out to be very efficient. 
These parametrizations allow for solving optimization problems within the fields 
of system identification and optimal control. 
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1. INTRODUCTION 

Lossless or stable allpass transfer functions play 
an important role in system theory mainly due 
to the Douglas-Shapiro-Shields factorization: any 
proper transfer function can be written as the 
product of a lossless function, which includes the 
dynamics of the system, and an unstable factor. 
In many problems in which a criterion must be 
optimized over a set of functions, the unstable 
factor of the optimum can be computed from 
the lossless one. This can be done in rational L 2 
approximation (Fulcheri and Olivi [1998]), system 
identification (Bruls et al. [1999]), multi-objective 
control (Scherer [2000]). These problems can thus 
be handled by optimization methods over the 
class of lossless functions of prescribed degree, 



or possibly a specified subclass. It is with such 
applications in mind that we will address the 
parametrization issue. 

An interesting and unusual approach of this opti- 
mization problems is to use the manifold structure 
of the class of lossless functions of fixed McMillan 
degree (Alpay et al. [1994]) and parameters com- 
ing from an atlas of charts. An atlas of charts 
attached with a manifold is a collection of lo- 
cal coordinate maps (the charts), whose domains 
cover the manifold and such that the changes of 
coordinates are smooth. Using such parameters 
allows to exactly describe the set on which an 
optimum is searched. This ensures that the op- 
timum will be stable and of the prescribed order. 
In practice, a search algorithm can be run through 



the manifold as a whole, using a local coordinate 
map to describe it locally and changing from one 
coordinate map to another when necessary. 

In the literature, atlases of charts have been de- 
rived both from the state-space approach using 
nice selections and from the functional approach 
using interpolation theory and Schur type al- 
gorithms. A connection between these two ap- 
proaches was found in the scalar (or SISO) case 
(Hanzon and Peeters [2000]) and generalized to 
the matrix case (Hanzon et al. [2004]). In this 
paper, an atlas is described in which balanced re- 
alizations can be computed from the Schur param- 
eters. The computation involves a product of uni- 
tary matrices and thus presents a nice numerical 
behavior. Moreover, for some particular choices of 
the interpolation points and directions, the bal- 
anced realizations possess a triangular structure 
which relates to nice selections (Peeters ct al. 
[2004]). 

The natural framework for these studies is that 
of complex functions. However, systems are often 
real- valued and their transfer functions T are real, 
that is, they satisfy the relation T{z) — T{z). 
Even if the complex case includes the real case 
by restriction, a specific treatment is actually rel- 
evant and was the initial motivation for this work 
which notably improves Marmorat et al. [2003]. 
In rational L 2 approximation for example, a real 
function may have a complex best approximant. 
This is the case for the function f(z) = 1/z 3 — 
1/z which admits three minima: a real one and 
two complex ones, which achieve the best relative 
error. 

In this paper, atlases are constructed in which 
lossless functions are represented by balanced 
realizations built recursively from interpolation 
data as in Hanzon et al. [2004]. But instead of 
the Nevanlinna-Pick interpolation problem used 
there we consider here the more general Nudelman 
interpolation problem. This very general frame- 
work allows to construct several atlases, includ- 
ing that of Hanzon et al. [2004], and to describe 
the subclass of real functions. For each partic- 
ular atlas presented in this work, we propose a 
simple method to find an "adapted chart" for a 
given lossless function. This last point, together 
with their nice numerical behavior, make these 
parametrizations an interesting tool for solving 
the optimization problems mentioned before. 



analytic in the open unit disk. To relate these two 
situations, we use the transformation F — > 
defined by 

F\z)=F*{l/z), F\z) = F{z)*. (1) 



Let 



J = 



I P 
-L 



A 2p x 2p rational matrix function Q(z) is called 
J-lossless (or conjugate J-inner) if, at every point 
of analyticity z of O(z) it satisfies 

<d(z)J<d(z)~* < J, |z| > 1, (2) 
<d(z)JQ(z)~* = J, \z\ = 1. (3) 

The simplest J-lossless functions are the constant 
J unitary matrices H satisfying H* JH = J. 

A p x p rational matrix function G(z) is called 
lossless or conjugate inner (resp. inner), if and 
only if 

G(z)G(z)* < I p , \z\ > 1 (resp. \z\ < 1), (4) 

with equality on the circle. The transfer function 
of a lossless system is a lossless function. A lossless 
function can have no pole on the unit circle and 
the identity G*(z)G(z) = I p for \z\ — 1, extends 
by analytic continuation to all points where both 
G(z) and G"(z) are analytic. Therefore, the func- 
tion G(z)' 1 agrees with G$(z) and is inner. 

We denote by the set of p x p lossless functions 
of McMillan degree n, by 7lL v n the subset of real 
functions and by U(p) the set of p x p constant 
unitary matrices. The McMillan degree will be 
denoted by deg. 

An important property of a lossless function is 
that if 

G(z) = C(zl n - A)~ 1 B + D, 

is a balanced realization (it always exists, see 
Genin ct al. [1983]), then the associated realization 
matrix 



R 



D C 
B A 



(5) 



2. PRELIMINARIES. 

This paper is concerned with finite dimensional, 
stable, discrete-time systems and their transfer 
functions which happen to be rational functions 
analytic outside the closed unit disk. Interpola- 
tion theory usually deals with functions that are 



is unitary. Lossless functions can thus be rep- 
resented by unitary realization matrices. Con- 
versely, if the realization matrix associated with 
a realization of order n of some p x p rational 
function G(z) is unitary, then G(z) is lossless of 
McMillan degree less or equal to n. For these 
questions, we refer the reader to Hanzon et al. 
[2004] and the bibliography therein. 



Along with a 2p x 2p rational function Q(z) block- 
partitioned as follows 



Q(z) 



&u(z) e 12 (z) 
&2i(z) e 22 (z) 



(6) 



with each block of size p x p, we associate the 
linear fractional transformation Tq which acts on 
p x p rational functions F{z) as follows: 



t q {f) = [e 11 F + e 12 ][e 2lJ F + e 2 



(7) 



For a composition of linear fractional transfor- 
mations, it holds that To o = Te*. Linear 
fractional transformations occur extensively in 
representation formulas for the solution of var- 
ious interpolation problems (Ball et al. [1990]). 
To adapt the results available in the literature 
for functions analytic in the disk to the case of 
functions analytic outside, we use the relation 



Q = T e (R) Q» = T Jl@Jl Ji = 



I, 



In particular, we have the following result, stated 
for inner functions for example in [Fulchcri and 
Olivi, 1998, Lemma 3]: 

Theorem 1. If Q(z) is a J- lossless matrix func- 
tion, then the map Tq sends every lossless func- 
tion to a lossless function. 



3. NUDELMAN INTERPOLATION FOR 
LOSSLESS FUNCTIONS 

The Nudelman interpolation problem is to find a 
pxp rational lossless function G(z) which satisfies 
an interpolation condition of the form 



^- J G t (z)U(zI s -Wy 1 dz = V, 



(8) 



where (U, W) is an observable pair and W is stable 
(U is p x S and W is <5 x S). Note that if W 
is a diagonal matrix, this problem reduces to a 
Nevanlinna-Pick problem. 

It is well-known that there exists a rational lossless 
function G(z) satisfying the interpolation con- 
dition (8) if and only if the solution P of the 
symmetric Stein equation 



P - W*PW = u*u - v*v 



(9) 



is positive definite [Ball et al., 1990, Th.18.5.2]. 
A triple (W, U, V) such that the solution P of (9) 
is positive definite, will be called an admissible 
Nudelman data set. A 2p x 2p J-lossless function 
can then be built from (W, U, V): 



@w,u,v(z) = 
[h P -{z-l)C{zI 5 



W) _1 P _1 (/ tf - w)-*c*f 



where C 



Theorem 2. Let (W, U, V) be an admissible Nudel- 
man data set, and let O = &w,u,vH, where 
®w,u,v is given by (10) and H is an arbitrary con- 
stant J-unitary matrix. For every lossless function 
F(z), the lossless function 



G = Tq (F) 



(11) 



satisfies (8) and deg G = deg F+S. Conversely, the 
set of all lossless solutions G(z) of (8) is given by 
(11) where F(z) is an arbitrary lossless function. 

Proof. This result is a particular case of [Ball 
et al., 1990, Th.18.5.2], which describes all the 
Schur functions which are solutions of a Nudelman 
interpolation problem. □ 

Let A and II be p x p unitary matrices. Then the 
following important relations are satisfied 



"A 


" 


@w,u,v 


'A* 


" 





n 





n* 



= 



W.AUMV 



W.AU.TIV 



{AF(z)W) = AT eWtUiV {F(z))U- 



(12) 
(13) 



4. BALANCED REALIZATIONS 



The aim of this section is to choose the arbitrary 
J-unitary factor H in Theorem 2, so that the 
linear fractional transformation G = Tq w u v h{G) 
yields a simple and powerful construction for 
balanced realizations as in (Hanzon et al. [2004]). 

Let U and V be (p + S) x (p + S) unitary matrices 
partitioned as follows: 



U 



M u a u 

K K u 



v = 



M v a v 



, (14) 



where K u and K v are 5 x 5, a u , a v , /3 U and /3 V are 
p x 8 and M u and M v are pxp, and put 



M = 



M u 

My 



a = 



OLv. 



(3 = 



.(15) 



Proposition 1. Let U and V be unitary matrices 
block-partitioned as in (14). Assume that K v z — 
K u is invertible. Given a p x p proper rational 
transfer function G(z) and a minimal realization 
G{z) = D + C(zl k - Ay 1 B, the formula 



D C 
B A 



U 

Jfc 



DOC 
Is 
BOA 



V* 






h 



(16) 



(10) 



in which D, D arc p x p, A is k x k, A is (S + k) x 
(5 + k), defines a mapping 

G(z) -»■ G(z) =D + C(zl s+k - A)- 1 ^. 



This mapping coincides with the linear fractional which gives 
transformation G — T<s> uv (G), associated with 
the 2p x 2p J-lossless function 



G{z) 1 = ($2iG + $22)(« , llG'+$is 



*w,v(«) = M + a{K v z - K u )- x ii*J 



I P 
zl v 



or equivalently G = T$ M V (G). It can be easily 
established that 



Proof. The case 5 = 1 has been studied in 
Hanzon et al. [2004]. It is easily verified that (16) 
defines a mapping since the function G(z) does not 
depend on the choice of the minimal realization 
(A, B, C, D) of G(z). 

We shall use a well-known formula for the inverse 
of a block matrix [Dym, 1989, sec. 0.2]. Assuming 
that the block d is invertible, the inverse of a block 
matrix is given by 



a b 
c d 



(18) 



(a*)- 1 -(o x )- 1 6d- 1 
-rf- 1 c(a x )" 1 [d- 1 +d- 1 c(a x )- 1 bd- 1 } 

where a x = a — bd~ 1 c is known as the Schur 
complement of a. In particular, if 



T(z) 



D C 

B {A- zl s+k ) 



= (1 — A z)a(K v z - K u y x {K v A - K u )~*a*, 

and thus $uy is ./-lossless. □ 

Remark. A state-space formula of the form (16) 
associated with some linear fractional transforma- 
tion has been used in Horiguchi [1999] to describe 
all the positive real functions which interpolate 
given input-output characteristics. 

Proposition 2. Let (W, U, V) be some admissible 
Nudelman data set. There exist unitary (p + 5) x 
(p + 5) matrices U and V and a 2p x 2p constant 
J-unitary matrix Hu y such that 







w,u,v Hu.y 



$>u,v 



(19) 



Proof. If (19) is satisfied, since &w,u,v(^) = hp, 
the matrix Huy must be given by 



Buy = M + a(K v - K u y l p*J. 



(20) 



then G(z)- 1 = [/,fl] T(z)- 1 



• By (16) 



T(Z)- 



U 

o h 



DOC 
Is 
BOA 



= r (*) 
r (z) = 



v* o 

h 

M U D 
'* U D - z t 
B 



' V* " 




'0 





. o I k _ 




_ zls+k _ 




M U C ' 




- zK v 









A - 


zlk_ 





The block matrix 
d = 



K u -zK v p* u C 
A-zIk 



is invertible and the Schur complement of M U D 
can be computed as 



M U D- [a u M U C] d- 1 
= $n(z)G(z) + $ 12 (z), 



B 



where $n(z) = M u - a u (K u - Kyz)^ 1 ^ and 
<&i2 (^) = +a u (K u — K v z)~ x fi*z are precisely the 
blocks of the function defined by (17). Still using 
(18) to compute r (z) _1 , we get 



G(z)- 1 = [M v a v 0] 



($ 11 G + $ 12 )- 1 
-rf- 1 c($ 11 G + $ 12 )- 



Moreover, the function $u,v cannot have a pole 
on the circle and can be rewritten 

^uy{z)H u 1 v = 

[hp -(z-l) a{K v z - K U )-\K V - K u )~*a*J] . 

The representation (10) of a J-lossless function 
being unique up to a similarity transformation 
(Ball et al. [1990]), there must exist a transforma- 
tion T such that P = T*T, K U K~ X = TWT^ 1 , 
ctK~ x = CT^ 1 . The matrix T is thus a square 
root of P. Since the matrix V must be unitary, we 
must have a^a v + K*K V = Is, or else 



[T~*V*VT~ 1 + Is)' 1 = K V K*. 



(21) 



The matrix T~*V*VT~ 1 + Is being positive defi- 
nite, this equation has solutions and K v being one 
of these, we can set 



in which 



a u = UK V 
K u = WK V , 
a v = VK V 



U = UT- 1 
W = TWT- 1 



(22) 



(23) 



V = VT~ 



Thcsc definitions imply that a^a u +K*K u = Is as 



required and the columns 



ot u 



and 



a v 



can 



be completed into unitary matrices U and V. The Consider the matrix 



columns 



B* 



and 



M v 



can be determined up 



to some right pxp unitary matrices. □ 



5. EXPLICIT FORMULAS FOR U AND V. 

An observable pair (U, W) such that W is stable 
is called output normal if it satisfies 



U*U + W*W = I S . 



(24) 



Note that two equivalent triples (W, U, V) and 
(TWT' 1 ,1/T- 1 ,VT~ X ) give the same interpola- 
tion condition (8), so that we can assume that the 
pair (U, W) in (8) is output normal. From now on, 
this normalization condition will be imposed to 
the admissible Nudclman data sets. Explicit for- 
mulas for U and V can then be given which ensure 
the smoothness of our parametrization. They have 
been used for implementation (see section 7.3). 

It has been proved that the second block columns 
of U and V are given by (22) and (23) in which T 
is a square root of P (the solution to (9)) and 
K v a solution to (21). We choose the uniquely 
determined Hermitian positive square roots T = 
P 1 ' 2 and K v = (I s + V*^)- 1 ' 2 . It remains to 
specify the completion of these block columns into 
unitary matrices. The matrix M v satisfies 

m v m: = i p - v{i s + v*^)- 1 ^*, 

and it is easily seen that 

i P - v{h + v *v)- x v* = (j p + vv*)- 1 

is positive definite. Thus we can choose 

My = (I p + VV*)- 1 ' 2 , 



I p 
T 



T~ 



X U 
TY W 



The problem is now to find a right factor of 
~* " 



the form 
matrix. Let 



★ K v 



which makes it into a unitary 



~ N L*~ 




X U ' 


* 


X U ' 


L K 




TY W_ 




TY W _ 



A classical method consists of writing a Cholesky 
factorization using the following well-known fac- 
torization of a (p + 6) x (p + S) block matrix 



N L* ' 




'l p VR- 1 ' 




' Z~ x " 




i P o 




L K 




h 









K 




K^L I S 






' z- 1 ' 2 




* 


z~ 


-1/2 " 








K- X ' 2 L K 1 ' 


2 






1/2 L R l/2 





where Z = (N — ^K^L)' 1 . By (18), Z is the 
left upper block of 
-i -l 



N L* 
L K 



I P 
T 



Ip o 

o p- 1 



Ip o 
T* 



and can be computed as 

Z = X*X + Y*P- 1 Y. 

The matrices L and K are given by 

L = U'*X + W*TY 
K = U*U + W*W = I s + V*V. 



(28) 



(29) 
(30) 



Note that the matrices Z and K are actually 
positive definite and that K^ 1 ' 2 = K v as desired. 
Thus, we can set 



U 



X U 
TY W 



Z 1 ' 2 
-K^LZ 1 ' 2 K- 1 ' 2 



(31) 



and thus /3y 
can set 



-V*{Ip + VV*) 1 ' 2 1 so that we This proves the following result: 



V 



{i p + vv*y 1/2 v(i s + v*vy 1/2 
-v*(i p + vv*)- 1/2 (Is + V*V)- 1/2 



.(25) 



The construction of a matrix Li is more involved 
since the matrix I p — a u o* u may fail to be positive 
definite. This is the case for example when W is 
the zero matrix, then a u = U and I p — UU* is not 
invertible. However, when V is zero, since W is 
stable, Is — W* is invertible and there is a simple 
way to construct a unitary matrix 



Proposition 3. Let (W, U, V) be some admissible 
Nudclman data set satisfying (24). Define the map 
r : (W, U, V) -> {U, V), where U and V arc the 
unitary matrices given by (31) and (25), where U, 
V, W are given by (23), K,L,Z by (30), (29), (28) 
and X, Y by (27), in all of which T is the positive 
square root of P, the solution to (9). 
Then, the J-losslcss function 



e 



w.u.y 



e 



w,u,vHuy, 



(32) 



U 



X u 
Y W 



X = I P -U(I S -W*)- X U*, 
Y = (Is -W)(I S -W*)- 1 !/*. 



(26) 



(27) 



Hu,v being given by (20), coincides with $>u,v- 

Let A, LT be p x p unitary matrices and S be a 
S x S unitary matrix. Noting that (AZA*) 1 / 2 = 
AZ 1 / 2 A*, it is easily verified that 



t{W, AU, UV) 
'A 
I s 



U 



A* 
I s 



n o 

Is 



V 



IT 

Is 



so that the J-lossless function @w,u,v also satis- 
fies (12) and (13). 

We also have that 
t(£*VF£,[/£,V/£) 



I P 

s* 



u 



Ip o 
E 



i P o 
s* 



V 



i P o 
s 



so that 



e 



s*w£.t/£.V£ 



e 



w,u,v- 



(33) 



Corollary 1. A unitary matrix realization R of 
G = Ta (G) can be computed from a uni- 
tary matrix realization R of G(z) by (16) in 
which U and V are the unitary matrices (U, V) = 
t(W, U, V) defined in Proposition 3. 



Remark. Instead of U and V, we may have chosen 



u = 





" 

°K 


w 


#1 




" 

o 2 


v = 


Ip 




" 

°1. 


V 


"#2 




" 

C-2 



in which Cq, O2, Hi and are unitary matrices: 

(1) the matrix 0\ corresponds to the choice 
of 0\P x l 2 instead of P 1 / 2 . This choice leave 
the linear fractional transformation unchanged 

= $w,v)j but changes the realization R = 
(A,B,C,D) into (OxAO^OxB, COl,D), a simi- 
lar one. 

(2) the choice of K v = (V*V + I p )- 1/2 2 instead 
of (V*V + Ip)- 1 ! 2 has no effect. 

(3) The unitary matrices H\ and H 2 correspond 
to another completion of the first columns of V 
and U. This choice changes the linear fractional 
transformation since 



uy 



uy 



Hi 
H 2 



and thus produces a non-similar realization. 

Remark. Note that if V = 0, since the pair 
{U, W) satisfies U*U + W*W = I s , we have that 
P = I s and T = I s too. Thus, U = U a defined 
by (26) and V = I p +s, so that the recursion (16) 
becomes 



D 


C 


B 


A 



' M U D 


U M U C~ 






w p* u c 


(34) 


B 


A 





6. CHARTS FROM A SCHUR ALGORITHM. 

Recall that a manifold is a topological space 
that looks locally like the "ordinary" Euclidean 
space R w : near every point of the space, we have 
a coordinate system or chart. The number N 
is the dimension of the manifold. It has been 
proved in [Alpay et al., 1994, Th.2.2] that C p n is a 
smooth manifold of dimension p 2 + 2np embedded 
in the Hardy space H PXp , for 1 < q < 00 

(TIC? is a smooth manifold of dimension + 
np). The topology on L p n is that induced by 
the L" norm on H pxp . In Alpay et al. [1994] 
atlases of charts have been constructed from a 
Schur algorithm associated with Nevanlinna-Pick 
interpolation. We generalize this construction to 
the case of Nudelman interpolation. 

Let a = ((I/i.Wi), (U 2 ,W 2 ),...,(Ui,Wi)) be a 
sequence of output normal pairs (see section 5), 
Wj is rij x rij , Uj is p x rij , and 



n. 



From a given lossless function G{z) of degree n, a 
sequence of lossless functions of decreasing degree 
Gi{z) = G(z), Gi-i(z), . . . can be constructed 
following a Schur algorithm: assume that Gj(z) 
has been constructed and put 



dz. 



If the solution P, to the symmetric Stein equation 



P j -W*P j W j = U*U j -V*V j 



is positive definite, then from Theorem 2, a loss- 
less function Gj-i(z) is defined by 



If Pj is not positive definite, the construction 
stops. 

A chart (D, cfi) of £P n is attached with a sequence 
a of output normal pairs and with a chart (W, tp) 
of U(p) as follows: 

A function G(z) G £ p belongs to the domain V 
of the chart if the Schur algorithm allows to con- 
struct a complete sequence of lossless functions, 

G(z) = Gi(z),G l - 1 (z)...,G , 

where Go is a constant lossless matrix in W C 
U(p). 

The local coordinate map <f> is defined by 

<j>:G(z)£V^(V 1 ,V 2 ,...,V l ,iP(G )), 



and the interpolation matrices Vj are called the 
Schur parameters of the function in the chart. 

Theorem 3. A family of charts (V, <ft) defines an 
atlas of L v n provided the union of their domains 
covers CL . 



for G(z). Such a chart presents a great interest 
from an optimization viewpoint. The optimization 
process starts, in an adapted chart, at the origin 
and thus far from the boundary where a change of 
chart is necessary. For each atlas, a simple method 
to find an adapted chart is given. 



Proof. The proof is analogous to the proof of 
Th.3.5 in Alpay et al. [1994]. The domain V of 
a chart is open and the map is a diffcomor- 
phism. This relies on the fact that Qw,u,v depends 
smoothly on the entries of V. □ 

Atlases for the quotient ££/U(p) are obtained 
using the properties (12) and (13). If G(z) has 
Schur parameters (Vj., V 2 , ■ . ■ , Vj) and constant 
unitary matrix Go in a given chart, and if 
II £ U(p), then G(z)U* has Schur parameters 
(nVi, IIV2, • • • , nV;) and constant unitary matrix 
Goll* in the same chart. The quotient can be 
performed within a chart by imposing the last 
constant lossless matrix Go in the Schur algorithm 
to be the identity matrix. 

In a chart, a balanced realization R of G(z) = 
-1 (Vi, V2, . . . , V'(Go)) can be computed from the 
parameters using the Schur sequence: let i? = 
Go, a realization Rj of Gj(z) is obtained from 
a realization Rj-i of Gj-i(z) by formula (16) 
in which (U,V) = (Uj,Vj) = T(Wj,Uj,Vj) (see 
Proposition 3). This process allows to select for 
each G(z) in the domain of the chart a unique 
balanced realization within the equivalence class, 
and then the map 

R^(V u V 2 ,...,iP(G )) 

is a canonical form. However, the domain of this 
canonical form is not easily characterized and it is 
in general difficult to decide if a given realization 
is in canonical form with respect to a chart. This 
can be done in some particular situations that will 
be studied in the following section. 



7. SOME PARTICULAR ATLASES. 

We describe three atlases which all present some 
interest from the optimization viewpoint. The 
first one is for complex functions and it involves 
only Schur steps in which the degree is increased 
by one. It allows for a search strategy of local 
minima by induction on the degree, which can 
be very helpful in some difficult optimization 
problems. The second one is the analog for real- 
valued functions. The third one involves only one 
Schur step and provides very simple and natural 
canonical forms. 

A chart in which all the Schur parameters Vj are 
zero matrices for G(z) is called an adapted chart 



7.1 An atlas for (complex lossless functions) 

Consider the charts associated with sequences of 
output normal pairs (m,wi), (u 2 , 102), ... , ( u n , w n ) 
in which the Wj's are complex numbers. In this 
case, the Nudelman interpolation condition (8) 
can be rewritten as a Nevanlinna-Pick interpola- 
tion condition 

G{l/wj)*Uj = Vj. 

This is the atlas described in Hanzon et al. [2004]. 
However, the normalization conditions differ. In 
Hanzon et al. [2004] the p-vectors Uj have norm 
one, while in this work, the pairs (uj,Wj) are 
output normal (24). 

Remark. Note that, when nj > 1 the normal- 
ization condition Uj unitary (a possible general- 



ization of 



= 1) cannot be chosen since the 



matrix U*Uj can be singular. 

In view of (34) , an adapted chart can be computed 
from a realization in Schur form. 

Lemma 1. Let G(z) £ C p n and let R = (A, B, C, D) 
be a balanced realization of G(z) in Schur form (A 
upper triangular). Let 



w a 
A 



B 



C=[uC], 



where w £ C, u, b £ <C P , and a £ (D™" 1 . Then, 
G = Tq (G) for some lossless function G(z). 
A realization R of G(z) can be computed by 
reverting (16). It is still in Schur form and given 
by R= (A,B,C,D), where 

G = G+(l-u;)-W 
D = D + {l-w)- 1 ub*. 

This process can be repeated. It provides a se- 
quence of output normal pairs (iij,Wj), the Wj's 
being the eigenvalues of A. In the corresponding 
chart the Schur parameters of G(z) are the zero 
p- vectors v n = . . . = v\ = 0. 



7.2 An atlas for TtZ v n (real lossless functions) 

To deal with real functions we consider the charts 
associated with sequences of output normal pairs 



(I/i.Wi), (U 2 ,W 2 ),...,(U n ,W n ) in which the 
Wj 's are either real numbers or real 2x2 matrices 
with complex conjugate eigenvalues and the Uj's 
are real matrices. The parameters (the matrices 
Vj) are then restricted to be real. 

As previously, an adapted chart for a given lossless 
function G(z) e Tl^ p n can be obtained from a 
realization in real Schur form. 



to (8) can be directly characterized in state space 
form and formula (16) recovered independently 
from the Schur algorithm. 

Proposition 4. Let G{z) = D + C(zI n -Ay 1 B be 
a balanced realization of G(z) G L v n . Let Q be the 
unique solution to the Stein equation 



Q - A*QW = C*U. 



(35) 



lemma 2. Let G{z) 6 TZ£p and let R = 
(A,B 7 C 7 D) be a balanced realization of G(z) in 
real Schur form 



A 



Wi * 

w x _ x 

••• 



• . * 
Wi 



where for j = 1, . . . , I, Wj is either a real number 
or a 2 x 2 block with complex conjugate eigenval- 
ues. Let 



Wi A* 
A 



C=[U l C], B 



B* 
B 



where Ui and B are pxni,ni being the size of W\. 



Then, G = Ta 



(G) for some lossless function 



G(z). A realization R of G(z) can be computed 
by reverting (16). It is in Schur form and given by 
R = (A,B,C,D), where 

C = C + {I ni -Wi)- 1 UiA* 
D = D + (I ni -Wi)- 1 U l B*. 

Repeating this process, we get a sequence of 
output normal pairs (Uj, Wj), the Wj's being the 
diagonal blocks of A, that index a chart in which 
the Schur parameters (which are now matrices of 
different sizes) of G(z) are all zero matrices. 

Remark. The Schur algorithm attached with 
G(z) in such an adapted chart yields a Potapov 
factorization for real lossless functions 

G(z)=B l (z)B l _ 1 (z)---B 1 (z), 



where Bj is the real- valued lossless function 



Bj{z) = 

I p -( z -l)Uj{zI nj -Wj)-\l nj 



W*) 



7.3 Lossless mutual encoding. 

For this atlas, we consider the charts associated 
with a single output normal pair (U, W) in which 
W is n x n and U is p x n. In this solution 



Then, the interpolation value V in (8) and the 
solution P to (9) are given by 



V = D*U + B*QW 1 
P = Q*Q. 

The unitary realization matrix R 
G(z) can be computed as 



D C 
B A 



(36) 
(37) 

of 



R = U 



Go 
/„ 



V*, Go e U(p), 



in which U and V are unitary matrices giyenas in 
Proposition 3 by (31) and (25), where U, V, W 
are given by (23), K,L,Z by (30), (29), (28) and 
A, Y by (27), but in which the square root T of P 
is now chosen to be Q. 

Proof. Since G*(z) = D* + B* {\l n - A*)^ 1 C* , 
the contour integral (8) can be computed as 



V =^~J GHz)U{zI n -W)- l dz 

T 

v =^l D ~ u {tr lw) ) 

I 00 \ / 00 

+b* z ^{zA*y c* u \J2 z ~ jwj 



dz 



0=0 



= d*u + b* ^J2(A*yc*uw j j w. 

Since Q is given by the convergent series Q = 
Y^L {A*yC*UW\ (36) is satisfied. 

Consider a unitary completion of the column 

^ , for example the matrix U given by (26). 

Here, 14$ and the unitary realization matrix R 
have the same size and formulas (35) and (36) 
can be rewritten in a matrix form 



R* 



I P 
Q 



* V 

* Q 



T. 



(38) 



The matrix R being unitary, we have that 



T*T=U* 



I P 
Q*Q 



U 



* V*V + Q*Q 



so that U*U + W*Q*QW = V*V + Q*Q, which 
proves (37). 



n(z) = x + u(zi n - w)- x y e cp. 



(39) 



We shall use the computations of section 4, with 
5 = n and T = Q instead of the Hermitian positive 
square root P 1 ! 2 . From (38) we get 



The canonical form associated with the pair 
(U, W) depends on this completion and is in fact 
attached with an element of ££/U(p). This ex- 
plains why this section was called lossless mutual 
encoding. 



and if IA is the unitary matrix given by (31), 
'* VK- 1 ' 2 ' 



R*U 



This matrix is unitary and its second block col- 
umn coincides with that of V given by (25). Thus 
it must be V up to a right unitary factor of the 



form 



h 



, for some unitary matrix Go 



Remark. Note that $^ v — <&uy (M, V defined in 
Proposition 3) and thus Go is the constant unitary 
matrix in the Schur algorithm G = T§ w u y (Go). 

The domain of a chart and of the associated 
canonical form can then be easily characterized. 



Proposition 5. A lossless function G{z), given 
by a balanced realization (A, B, G, D) can be 
parametrized in the chart defined by the pair 
(U, W) if and only if the solution Q to the Stein 
equation (35) is positive definite. A realization R 
is in canonical form with respect to this chart if 
and only if the solution Q to (35) is P 1 / 2 , P being 
a solution of (9). 

The invertibility of the matrix Q is a good mea- 
sure of the quality of the chart, the best choice 
being Q — I n . This choice provides an adapted 
chart. 



8. APPLICATION TO SYSTEM 
IDENTIFICATION AND CONTROL. 

The first application that we consider is the identi- 
fication of hyperfrequency filters, made of coupled 
resonant cavities, that are used in telecommunica- 
tion satellites for channel multiplexing. The prob- 
lem is to recover the transfer function of the fil- 
ter from frequency data. These data are estimate 
values of the transfer function at pure imaginary 
points obtained from the steady-state outputs of 
the filter to harmonic inputs. A first stage, far 
from being trivial, consists in computing a sta- 
ble matrix transfer function of high order which 
agrees with the data. It is achieved by the software 
Hyperion, also developed at INRIA (Baratchart 
et al. [1998]). Then a rational L 2 approximation 
stage is performed by the software RARL2 1 , in 
which the atlas of section 7.1 is used. Transfer 
functions of these filters are complex functions 
since a particular transformation has been used to 
simplify the model. In Figure 1, a 8th order model 
of a MIMO 2x2 hyperfrequency filter is shown, 
obtained from 800 pointwise data. A longstanding 



rr=0.63% err.sup =1.0 




ditf-12-21 =4.57% 



Proposition 6. The chart associated with the out- 
put normal pair (G, A) in a balanced realization 
(A, B, G, D) of G(z) is an adapted chart for G(z). 



r=0.04% err=1.22% err.sup =1.90% diff-12-21 =4.57% 



r=0.33% err=0.69%, err.sup =1 .24% 



Proof. In this case, Q = I, 
V = 0. 



n , so that P = I n and 
□ 



Remark. Let (£>», : G(z) ->■ (V,G ) be 
the chart associated with the output normal pair 
(U, W). Then the chart associated with the pair 
(UE,T.*WT.) has same domain V and by (33) 
coordinate map <j)' : G(z) — > {VH, Go). In an atlas, 
these two charts play the same role. 

Remark. Equivalence classes of output normal 
pairs are in bijection with lossless functions in 
[Alpay et al., 1994, Cor.2.1]. The unitary 

completion Uq of the matrix 1 



a lossless function 



IF 



in (26) defines 




Fig. 1. CNES 2x2 hyperfrequency filter: data and 
approximant at order 8 (Bode diagram). 

cooperation with the space agency CNES resulted 
in a dedicated software PRESTO-HF that wraps 
both HYPERION and RARL2 into a package 
which is now fully integrated in the design and 
tuning process. 



1 The software RARL2 is described in Marmorat et al. 
[2002] and available at the web page 

http: www-sop. inria. fr/apics/RARL2/rarl2-eng. html . 



This application and more generally filter de- 
sign, raises interesting new parametrization is- 
sues. The physical laws of energy conservation 
and reciprocity introduce subclasses of transfer 
functions which play an important role in this 
domain. These include J-inner, Schur (or contrac- 
tive), positive real and symmetric functions. In 
another connection, systems having a particular 
state space form must be handled to account for 
some physical properties, like for example the 
coupling geometry of a filter. We think that Schur 
analysis could help us to describe such subclasses 
and to pave the bridge between the frequency 
domain (where specifications are made) and the 
state-space domain (where the design parameters 
live) . As a first step in this way, a Schur algorithm 
for symmetric lossless functions, based on a two- 
sided Nudelman interpolation condition, has been 
presented in Olivi et al. [2005]. 

The atlas of section 7.3 has been recently imple- 
mented in RARL2, the rational L 2 approximation 
software. Its effectiveness has been demonstrated 
on random systems and classical examples from 
the literature. A promising field of application, in 
which functions are real- valued, is multi-objective 
control. In Scherer [2000], revisited in a chain- 
scattering perspective in Drai et al. [2005], it is 
shown that if the pair (Cq,Aq) of the Youla pa- 
rameter Q(z) = Dq + Cq(zI — Aq)~ 1 Bq is fixed, 
then the search over the parameters (Bq, Dq) can 
be reduced to an efficiently solvable LMI problem. 
Limiting the search of the parameter Q(z) to the 
FIR form 

Q(z) = Qo + Qi- + ... + Q P 4 

z zP 

provides solutions to the multi-objective control 
problem. However, this is also the main limitation 
of the approach as high order expansions might be 
necessary, due notably to the fact that the poles 
structure is fixed through the pair (Cq, Aq). Such 
a drawback could be avoided if the search was 
performed over all the parameters Q(z) of fixed 
McMillan degree. This could be done using the 
atlas of section 7.3 to parametrize the correspond- 
ing pairs (Cq,Aq). This work is currently under 
investigation and will be reported later. 
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